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Abstract 

The objective of this research is to improve and implement the filtered mass 
density function (FDF) methodology for large eddy simulation (LES) of high-speed 
reacting turbulent flows. We have just completed Year 1 of this research. This is the 
Final Report on our activities during the period: January 1, 2003 to December 31, 
2003. 2002. 

In the efforts during the past year, LES is conducted of the Sandia Flame D 
[1,2], which is a turbulent piloted nonpremixed methane jet flame. The subgrid 
scale (SGS) closure is based on the scalar filtered mass density function (SFMDF) 
methodology [3]. The SFMDF is basically the mass weighted probability density 
function (PDF) of the SGS scalar quantities [4]. For this flame (which exhibits 
little local extinction), a simple flamelet model is used to relate the instantaneous 
composition to the mixture fraction. The modelled SFMDF transport equation is 
solved by a hybrid finite-difference/Monte Carlo scheme. 

Dr. J. Philip Drummond (Hypersonic Propulsion Branch, NASA LaRC, Mail Stop 
197, Tel: 757-864-2298) is the Technical Monitor of this Grant. 



1 Introduction 


There have been significant recent developments of subgrid scale (SGS) clo- 
sures for large eddy simulation (LES) of turbulent reacting flows. Several re- 
cent reviews are available [4-10]. One such closure is via the filtered density 
function (FDF) methodology, first introduced by Pope [11]. This is the coun- 
terpart of the probability density function (PDF) method which has proven 
quite effective in Reynolds averaged simulations (RAS) [4,12]. This success 
is due to the inherent property of the PDF of providing complete statistical 
information about the variables. Due to this property, the FDF offers the abil- 
ity to treat finite-rate chemistry and the turbulence-chemistry interactions. In 
comparison to RAS, the LES/FDF methodology provides a more detailed and 
reliable prediction of turbulent reacting flows and is better suited to account 
for the large scale unsteady phenomena which are prevalent in combustion 
devices. The scalar FDF (SFDF) is considered by many investigators [13-17]. 
Its extension to account for variable density flows is via the scalar filtered mass 
density function (SFMDF) [3,18]. The velocity FDF (VFDF) is introduced by 
Gicquel et al. [19] and the velocity-scalar FDF (VSFDF) by Sheikhi et al. [20]. 
A recent review on the state-of- progress on LES/FDF is provided by Givi [10]. 
The outcome of these investigations has been encouraging, warranting further 
extension and application of this methodology in turbulent combustion. 

In this work, we employ the FDF method for prediction of the piloted jet 
flame studied in the experiments of the Combustion Research Facility at the 
Sandia National Laboratories [21,22]. This flame has been the subject of broad 
investigations by other computational/modelling methodologies [2,23-26]. In 
the experiments, three turbulent flames are considered: Flames D, E and F. 
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The geometrical configuration in these flames is the same, but the jet inlet 
velocity is varied. In Flame D, the fuel jet velocity is the lowest and the flame 
is close to equilibrium. The jet velocity increases from flames D to E to F, with 
noticeable non-equilibrium effects in the latter two. Flame D is considered in 
this work. The objective is to assess the predictive capability of the LES/FDF 
methodology in capturing the flow field and scalar mixing. This is a necessary 
step before consideration of the non-equilibrium flames (E and F). 


2 Formulation 

Sandia Flame D consists of a main jet with a mixture of 25% methane and 
75% air by volume. The nozzle is placed in a coflow of air and the flame 
is stabilized by a substantial pilot. The Reynolds number for the main jet 
is Re = 22400 based on the nozzle diameter D — 7.2 mm and the bulk jet 
velocity 49.6 m/sec. The methane-air reaction mechanism, as occurs in this 
flame, is taken into account via the “flamelet” model. This model considers 
a laminar, one-dimensional counterflow (opposed jet) flame configuration [6]. 
The detailed kinetics mechanism of the Gas Research Institute (GRI2.11) [27] 
is employed to describe combustion. At low strain rates, the flame is close to 
equilibrium. Thus, all of the thermo-chemical variables are determined by the 
“mixture fraction.” 

Formulation is based on the compressible form of the continuity, Navier- Stokes, 
energy (enthalpy) and mixture fraction conservation equations in a low Mach 
number flow [28]. These equations govern the space (x = x i} i = 1, 2, 3) and 
time ( t ) variations of the fluid density p(x, t), the velocity vector u — u*(x, f), 
the pressure p(x, t), the specific enthalpy h(x,t) and the mixture fraction 
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£(x, t). We employ Fourier’s law of heat conduction, Fick’s law of diffusion 
and we assume unity Lewis number. The molecular viscosity n increases with 
temperature (T) to the power of 0.7. The magnitude of the molecular Schmidt 
(and Prandtl) number is Sc = 0.75. 

For LES, we consider a spatially and temporally invariant positive filter £/(x' — 
x) of width Ag. The filtered value of the variable Q(x, t) is denoted by 
(Q(x,t)) e and its mass- weighted filtered value is denoted by (Q(x, t)) L = 
( pQ)i / ( p)e ■ The SFMDF of the mixture fraction is denoted by x ,t) 

where ip denotes the (probability) composition domain of the mixture frac- 
tion. The SFMDF accounts for SGS statistics of only the scalar field. Thus, 
the closure of the SGS hydrodynamics must be provided by other means. 
Here, we employ a well-established gradient diffusion model, in which the 
SGS dynamic viscosity u t is modelled by the MKEV (modified kinetic energy 
viscosity) [3,29] model. In the transport equation governing the SFMDF, the 
effects of SGS convection are also modelled by an analogous gradient diffu- 
sion model. The SGS mixing is closed via the least square mean estimation 
(LMSE) model [30,31]. These are described in detail in our previous papers 
on SFDF [3,15]. The final, modelled SFMDF transport equation reads: 


dt dxi 


_d_ 

dxi 


(7 + 7t) 


{p),y 

dxi 




(i) 


Here, 7 = ^ is the molecular diffusivity, is the SGS diffusivity 

and Sc t = 0.75 is the SGS Schmidt number. With the closure as such, the 
modelled scalar flux is consistent with that in most other (non-FDF) methods 
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[32,33]: 


The term Cl m is the SGS mixing frequency (as required in LMSE) and is 
modelled by [15,3] = Cq ( 7 + 7 t)/((p) e A|), with Cq = 8. 

The relation between the thermo-chemical variables (denoted by the array 
</>(x, f)) and the mixture fraction is provided by the flamelet library. In the 
context of the opposed jet flame, the library provides: 4> = </>(£, a) where 
a denotes the strain rate. For a fixed value of a, the SGS statistics of the 
thermo-chemical variables are determined from the SFMDF: 

(0(x,f)) L = *(0^(6x,t)«. (3) 


A hybrid finite-difference (FD) / Monte Carlo (MC) scheme is employed to 
solve the coupled set of the SFMDF equation and the filtered hydrodynamics 
equations. In this scheme, the domain is discretized by FD grid points and 
the SFMDF is represented via an ensemble of MC particles [34]. All of the 
hydrodynamic variables are determined on the FD points. A fourth order 
compact scheme [35,36] is used for FD discretization. Transport of the MC 
particles and the change in their properties are described by a set of stochastic 
differential equations (SDEs) similar to those describing a diffusion process 
[37]. The MC particles undergo motion in physical space by convection due to 
the filtered flow velocity and diffusion due to molecular and SGS diffusivities. 
The compositional value of each particle is changed due to SGS mixing. These 
are described by the SDEs: 


dxf = 


( u *)l + 


1 5(7 + 7t) 


(p)e dxi 


dt + 1/2(7 + 7t)/ (p)e dwu 


de = -a» [£ + ~ <0 l] <**> 
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(4) 

(5) 



where xf and £ + denote Lagrangian position and the composition (mixture 
fraction), respectively. The term w denotes the Wiener process [38]. The 
Fokker-Planck equation corresponding to Eqs. (4)-(5) is equivalent to Eq. (1). 
Thus, the solution of these SDEs represent the FMDF in the probabilistic 
sense. 


To understand the operational procedure, the elements of the computation 
are shown in Fig. 1. For simplicity, a two-dimensional domain is shown with 
fixed grid spacing of size of A. The MC particles are distributed randomly 
and are free to move within the domain as governed by Eq. (4). This trans- 
port is Lagrangian, thus the solution is free of the constraints associated with 
typical simulation of convection on fixed grid points. Statistical information, 
e.g. filtered values, at any point is obtained by considering an ensemble of Ne 
computational particles residing within an ensemble domain of side length A e 
centered around the points. For reliable statistics with minimal numerical dis- 
persion, it is desired to minimize the size of ensemble domain and maximize 
the number of the MC particles [12] . In this way, the ensemble statistics would 
tend to the desired filtered values 


<«» ■ w B £ 


Ne — >oo 
A ^— >0 


W)z. 


(6) 


where Q ^ denotes the information carried by n th MC particle pertaining to 
transport variable Q. Transfer of information from the grid points to the MC 
particles is accomplished via interpolation. The transfer of information from 
the particles to the grid points is accomplished via ensemble averaging. 


In the hybrid scheme, some of the filtered quantities are obtained by MC, some 
by FD and some by both. That is, there is a “redundancy” in determination 
of some of the quantities. This redundancy is very useful in monitoring the 


5 



consistency of the simulated results [39]. In general, all of the equations for 
the filtered quantities can be solved by FD, where all of the unclosed terms 
are evaluated by MC. This process can be done at any filtered moment level 
[19]. Here, the filtered values of the mixture fraction and the temperature are 
used to establish consistency. Some results in this regard are presented in the 
next section. In addition, the computational accuracy of the methodology is 
established by examining both the statistical and the dispersion errors [40]. 
In doing so, the correlation of the fluid density with distribution of the MC 
particles, the size of the ensemble domain, the number of particles within this 
domain and the global distribution of the particles are monitored in a manner 
similar to that reported in our previous work [3,15,19,20]. 

3 Results 

In the simulations, all of the flow variables at the inflow are set the same as 
those in the experiments, including the inlet profiles of the velocity and the 
mixture fraction. The flow is excited by superimposing oscillating axisymmet- 
ric and helical perturbations at the inflow. The procedure is similar to that 
of Danaila and Boersma [41], but is modified to match the experimental data 
in the near field. Simulations are conducted on a three-dimensional Cartesian 
mesh with uniform spacings in each of the three directions. The computational 
domain spans a region of 18Z? x 10D x 10 D in streamwise, and the two lateral 
directions. The number of grid points are 91 x 101 x 101, respectively. The 
filter size is set equal to A g = 2(AxAyAz) ( - 1 ^ where Ax, Ay and Az denote 
the grid spacings in streamwise ( x ) and the two lateral ( y — z ) directions, 
respectively. The size of the ensemble domain for evaluations of the filtered 
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values from the MC solver is equal to the filter size. There are approximately 
48 MC particles at each grid point for these evaluations. In total, there are 
about 3.4 million MC particles within the domain at all times. 

The flamelet table at strain rate of a = 100 1/s is used to relate the thermo- 
chemical variables to the mixture fraction. First, the consistency and accuracy 
of the computations are assessed. Next, the overall predictive capability of 
LES/SFMDF is demonstrated by comparing the flow statistics with the Sandia 
data. These statistics are obtained by long-time averaging of the filtered field. 
The notations Q and RMS(Q) denote, respectively, the time-averaged mean 
and root mean square values of the variable Q. 

The simplest consistency check is via flow visualization. For example, Fig. 2 
shows the instantaneous contours of the filtered temperature field as obtained 
by the FD and the MC methods. The central jet lies in the middle along the 
axial coordinate, surrounded by a pilot where the temperature is the high- 
est and encircled by the air coflow. The region close to inlet is dominated by 
the molecular diffusion and the jet exhibits a laminar-like behavior. Farther 
downstream, the growth of perturbations is manifested by the formation of 
large scale coherent vortices. The upstream feedback from the vortices created 
initially is sufficient to trigger self-sustaining vortex rollup and subsequent 
pairing and coalescence of neighboring vortices. Due to presence of helical in- 
stabilities, the layer is asymmetric. The similarity of the results in the two 
figures is observed at all other times and is also observed for the mixture frac- 
tion (not shown) . This consistency is further assessed via establishing identical 
statistics of the redundant quantities as generated by the two methods. 

The capability of the method in predicting the hydrodynamics field is demon- 
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strated by examining some of the (reported) flow statistics. The centerline 
mean and RMS values of the axial velocity are compared to experimental 
data in Fig. 3. This figure indicates that the flow is adequately excited and 
the predicted results are in good agreement with data. The statistics of the 
thermo-chemical variables are also compared with corresponding data. The 
measurements are exhaustive; here, only some sample results are shown. The 
radial (r = \/z^~+ y 2 ) distribution of the mixture fraction is shown to compare 
well with data (Fig. 4a). Similar agreement is observed at other streamwise 
locations. The mean temperature values in Fig. 4b indicate a slight over- 
prediction on the rich side. This is due to premixing of methane with air as 
also indicated previously [25] . The “resolved” RMS values of the mixture frac- 
tion and temperature are in good agreement with data. However, the “total” 
RMS values, including the contributions of both the resolved and the SGS 
fields, are higher than values reported experimentally. The contribution of the 
SGS to the total scalar energy is about 30% which is expected in LES. The 
higher values of the total RMS values as predicted by LES/FD are not due 
to MC numerical dispersion because the FD results do indeed yield the same 
values. The level of SGS variance can be decreased by increasing the magni- 
tude of Cq. However, this would not alter the total RMS values significantly. 
It is possible that some contribution of this variance is not included in the 
measurements due to finite probe size. Higher resolution measurements would 
determine the allocations of scalar variance to the resolved and SGS fields. 

The statistics of the mass fractions (denoted by Y) of several of the species are 
compared with data in Fig. 5. Similar to the temperature results on the rich 
side, the reactants’ mass fractions are under-predicted slightly, while the prod- 
ucts are over-predicted. The mean values of the mass fractions of the major 
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and the minor species compare well with experimental data. All of the results 
indicate the adequacy of the flamelet table in relating the thermo-chemical 
variables to the mixture fraction, and also the good predictive capability of 
the SFMDF for this flame. The level of agreement of the RMS values of the 
mass fraction is the same as that presented in Fig. 4. 

The computational costs associated with LES/FDF depends, obviously, on 
the parameters of the simulations. For the case reported here, the simulations 
required about 110 hours of CPU time on a SUN Fire 4800 with 6 processors. 
This includes the times required for consistency tests and ensemble averaging 
of data. The computational time for LES without including SGS effects is 
about 10-12 times less. However such simulations yield erroneous predictions 
and in many cases lead to numerical instabilities. For further comparative 
assessment of the computational requirements of the FDF in comparison to 
non-FDF methods, we refer to our previous work[3, 15, 19,20]. 


4 Summary and Concluding Remarks 

The filtered density function (FDF) methodology [11] is now at a stage that it 
can be used for prediction of complex turbulent reacting flows. This is demon- 
strated in this work by utilizing the simplest form of the FDF for large eddy 
simulation (LES) of a piloted, nonpremixed, turbulent, methane jet flame 
(Sandia Flame D). For this near-equilibrium flame, the thermo-chemical vari- 
ables are related to the mixture fraction. This is done by construction of a 
flamelet library (in a counter-flow jet flame) in which the chemical reaction 
is modelled by detailed kinetics [27]. It is useful to note that the approach 
here is fundamentally different from those followed in previous flamelet based 
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SGS models. In most previous contributions [25,42-45], the FDF of the mix- 
ture fraction is assumed ( e.g . beta or other distribution). Here, a modelled 
transport equation for the mass weighted FDF of the mixture fraction [3] is 
considered. This equation is solved by a hybrid finite-difference / Monte Carlo 
method. After establishing the consistency and accuracy of the hybrid solver, 
the predictive capability of the overall scheme is assessed by comparison with 
experimental data. For these comparisons, the ensemble (long time averaged) 
values of the thermo-chemical variables are considered. It is shown that all 
of the mean quantities are, generally, predicted well. The resolved RMS of 
the variables also compare well with data. However, when the contribution 
of the subgrid scale (SGS) quantities are included, the experimental data are 
over-predicted. 

There are two ways by which this work can be continued. First, is extension 
of LES/SFMDF for prediction of flames which experience extinction (such as 
Sandia Flames E and F) and/or re- ignition. This would provide a more defini- 
tive assessment of the predictive capabilities of the FDF. Such simulations 
require consideration of finite-rate chemistry, as demonstrated in RAS/PDF 
simulations of Sandia flames [23,24]. Presently, it is not computationally fea- 
sible to implement detailed kinetics in such simulations. But implementation 
of reduced kinetics schemes using in situ adaptive tabulation, such as those 
used in RAS/PDF [23,24] is within reach. Second, it is desirable to implement 
the LES/VSFDF [20] for prediction of this (or other complex) flame(s). In 
VSFDF, the SGS convection appears in a closed form, thus the assumption of 
gradient diffusion can be relaxed. This is also feasible within the near future. 
Accomplishments of both of these tasks can be expedited by further reduction 
of the computational requirements of the FDF. The present work establishes 
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the capability of LES/FDF for accurate prediction of complex flames, war- 
ranting its further use for modelling of even more complex turbulent reacting 
flows. 
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Fig. 1. Elements of computation as used in a typical LES/FDF. Solid squares denote 
the FD points and the circles denote the MC particles. Also shown are three different 
ensemble domains. 
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Fig. 2. Instantaneous contours of the filtered values of the temperature ( K ) as 
obtained by FD and MC simulations. 
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x/D 

Fig. 3. The mean and RMS values of the axial velocity [m/s] at the centerline. Ucl 
denotes the mean axial velocity at the centerline, the symbols denote experimental 
data and the lines denote the predictions. 
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Fig. 4. Radial distributions of the mean and RMS values at x/D = 15. The symbols 
denote experimental data. The thick lines denote the mean values, the thin solid 
lines denote the resolved RMS values and the thin dashed lines denote total RMS 
values, (a) Mixture fraction, (b) temperature. 
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Fig. 5. Radial distribution of the mean values of the mass fractions at x/D = 15. 
The symbols denote experimental data, (a) Solid line: CH^\ dashed line: 02? (b) 
solid line: CO 2 ; dashed line: CO. 
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